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Abstract 

The paper is devoted to a constitutive solution, limit load analysis and Newton-like methods 
in elastoplastic problems containing the Mohr-Coulomb yield criterion. Within the constitu¬ 
tive problem, we introduce a self-contained derivation of the implicit return-mapping solution 
scheme using a recent subdifferential-based treatment. Unlike conventional techniques based 
on Koiter’s rules, the presented scheme a priori detects a position of the unknown stress tensor 
on the yield surface even if the constitutive solution cannot be found in closed form. This 
fact eliminates blind guesswork from the scheme, enables to analyze properties of the constitu¬ 
tive operator, and simplihes construction of the consistent tangent operator which is important 
for the semismooth Newton method applied on the incremental boundary value elastoplastic 
problem. The incremental problem in Mohr-Coulomb plasticity is combined with the limit load 
analysis. Beside a conventional direct method of the incremental limit analysis, a recent indirect 
one is introduced and its advantages are described. The paper contains 2D and 3D numerical 
experiments on slope stability with publicly available Matlab implementations. 

Keywords: infinitesimal plasticity, Mohr-Coulomb yield surface, implicit return-mapping scheme, 
consistent tangent operator, semismooth Newton method, incremental limit analysis, slope stability 

1 Introduction 

This paper is a continuation of [1] which was devoted to a solution of elastoplastic constitutive 
problems using a subdifferential formulation of the plastic flow rule. It leads to simpler and more 
correct implicit constitutive solution schemes. While a broad class of elastoplastic models containing 
1 or 2 singular points (apices) on the yield surface was considered in [1], the aim of this paper 
is to approach the subdifferential-based treatment to models that are usually formulated in terms 
of principal stresses. For example, the principal stresses are used in models containing the Mohr- 
Coulomb, the Tresca, the Rankine, the Hoek-Brown or the unihed strength yield criteria lassiE]. 
Such criteria have a multisurface representation leading to a relatively complex structure of singular 
points. 

Due to technical complexity of implicit solution schemes for these models, we focus only on a 
particular but representative yield criterion: the Mohr-Coulomb one. This criterion is broadly ex¬ 
ploited in soil and rock mechanics and its surface is a hexagonal pyramid aligned with the hydrostatic 
axis (see, e.g., 0)- We consider the Mohr-Coulomb model introduced in [21 Section 8] which can 
optionally contain the nonassociative flow rule and the nonlinear isotropic hardening. The nonasso- 
ciative flow rule enables to catch the dilatant behavior of a material. Further, due to the presence of 
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the nonlinear hardening, one cannot hnd the implicit constitutive solution in closed form, and thus 
the problem remains challenging. As in [2], we let a hardening function in an abstract form. For a 
particular example of the nonlinear hardening in soil mechanics, we refer, e.g., [6]. 

In literature, there are many various concepts of the constitutive solution schemes for models 
containing yield criteria written in terms of the principal stresses. For their detailed overview and 
historical development, we refer the recent papers [3] and [7], respectively. It is worth mentioning 
that the solution schemes mainly depend on a formulation of the plastic flow rule, its discretization 
and other eventual approximations. 

In engineering practice, the plastic flow rule is usually formulated using the so-called Koiter rule 
introduced in ^ for associative models with multisurface yield criteria. Consequently, this rule was 
also extended for nonassociative models, see, e.g., [S]. It consists of several formulas that depend 
on a position of the unknown stress tensor cr on the yield surface. The formulas have a different 
number of plastic multipliers. Within the Mohr-Coulomb pyramid, one plastic multiplier is used for 
smooth portions, two multipliers at edge points, and six multipliers at the apex. For each Koiter’s 
formula, a different solution scheme is introduced. However, only one of which usually gives the 
correct stress tensor. Moreover, the handling with different numbers of plastic multipliers is not 
suitable for analysing the stress-strain operator even if the solution can be found in closed form. 
If an elastoplastic model contains a convex plastic potential as the Mohr-Coulomb one then it is 
possible to replace the Koiter rule with a sub differential of the potential (see, e.g., m)- Such a 
formulation is independent of the unknown stress position, contains just one plastic multiplier, and 
thus it is more convenient for mathematical analysis of the constitutive operators. In [T], it was 
shown that this formulation is also convenient for a solution of some constitutive problems. Further, 
in some special cases, the constitutive problem can be also dehned using the principle of maximum 
plastic dissipation [21 |T0] or by the theory of bipotentials m and solved by techniques based on 
mathematical programming. 

We focus on the (fully) implicit Euler discretization of the flow rule, which is frequently used 
in elastoplasticity. Beside other Euler-type methods (see, e.g., mm), the cutting plane methods 
are also popular. We refer, e.g., na for the literature survey and recent development of these 
methods. When the constitutive problems are discretized by the implicit Euler methods, the solution 
is searched by the elastic predictor - plastic correction method. Within the plastic correction, the so- 
called (implicit) return-mapping scheme is constructed. It is worth mentioning that plastic correction 
problems can be reduced to problems formulated only in terms of the principal stresses [iaEii2|. 

In order to simplify the solution schemes for nonsmooth yield criteria, many various approximative 
techniques have been suggested. These techniques are based on local or global smoothing of yield 
surfaces or plastic potentials. For literature survey, we refer [3l Section 1.2] or [T5l[T6l[6]. However, 
such an approach is out of the scope of this paper. 

The constitutive problem is an essential part of the overall initial boundary value elastoplastic 
problem. Its time discretization leads to the incremental boundary value problem which is mostly 
solved by nonsmooth variants of the Newton method [T71 ITSl fl^ l2n[ l2T] in each time step. Then, it 
is useful to construct the so-called consistent tangent operator representing a generalized derivative 
of the discretized constitutive stress-strain operator. We use the framework based on the eigenpro- 
jections of symmetric second order tensors, see, e.g., [221 m. A similar approach is also used in 
the recent book [16] with slightly different terminology like the spectral directions or the spin of a 
tensor. Another approach is introduced, e.g., in [23l [111 E] where the consistent tangent operator is 
determined by the tangent operator representing the relation between the stress and strain rates. 

Further, this paper is devoted to the limit load problem which is frequently combined with the 
Mohr-Coulomb model. It is an additional problem to the elastoplastic one where the load history 
is not fully prescribed. It is only given a hxed external force that is multiplied by a scalar load 
parameter whose limit value is unknown. It is well known that the investigated body collapses 
when this critical value is exceeded. Therefore, this value is an important safety parameter and 
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beyond it no solution exists. Strip-footing collapse or slope stability are traditional applications on 
this problematic (see, e.g. mw- The simplest computational technique is based on the so-called 
incremental limit analysis where the load parameter is enlarged up to its limit value. Then, the 
boundary-value elastoplastic problem is solved for investigated values of this parameter. Beside the 
conventional direct method of the incremental limit analysis, we also introduce the indirect method 
and describe its advantages based on recent expertise introduced in [23 123 [23,12B] • 

The rest of the paper is organized as follows. In Section an auxilliary framework related to the 
subdifferential of an eigenvalue function and derivatives of eigenprojections is introduced. In Section 
[3 the Mohr-Coulomb constitutive initial value problem is formulated using the subdifferential of 
the plastic potential and discretized by the implicit Euler method. In Section the existence and 
uniqueness of a solution to the discretized problem is proven and the improved solution scheme is 
derived. In Section the stress-strain and the consistent tangent operators are constructed. In 
Section the direct and indirect methods of the incremental limit analysis are introduced. Both 
methods are combined with the semismooth Newton method. In Section [3 2D and 3D numerical 
experiments related to slope stability are introduced. In Section some concluding remarks are 
mentioned. The paper also contains Appendix with some useful auxilliary results. In Appendix 
A, the solution scheme is simplihed under the plane strain assumptions. In Appendix B, algebraic 
representations for second and fourth order tensors within the 3D and plane strain problems are 
derived. 

In this paper, second order tensors, matrices, and vectors are denoted by bold letters. Further, 
the fourth order tensors are denoted by capital blackboard letters, e.g., De or I. The symbol ( 8 ) means 
the tensor product |2]. We also use the following notation: M_|_ := { 2 ; G M; z > 0} and for 
the space of symmetric, second order tensors. The standard scalar product in and the biscalar 
product in are denoted as ■ and :, respectively. 

2 Subdifferentials and derivatives of eigenvalue functions 

In this section, we introduce an auxilliary framework that will be crucial for an efficient construction 
of the constitutive and consistent tangent operators in Mohr-Coulomb plasticity. Let 
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'n = '^ ® e*, T]i>r]2> 773, ( 2 . 1 ) 

i=l 

be the spectral decomposition of a tensor rj G Here, rji G M, e* G i = 1,2,3, denote the 

eigenvalues, and the eigenvectors of 77 , respectively. The eigenvalues rji, 772 , Vs can be computed using 
the Haigh-Westargaard coordinates (see, e.g., [21 Appendix A]), and they are uniquely determined 
with respect to the prescribed ordering. Let ui,U2,cu3 denote the corresponding eigenvalue functions, 
i.e. Vi •= ^ = 1,2,3. Further, we dehne the following set of admissible eigenvectors of r/; 

1 /( 77 ) = {( 61 , 62 , 63 ) G X X I 6 i ■ ej = dij] r]ei = r]iei, ij = 1,2,3; Vi > V 2 > Vs}- 

2.1 Subdifferential of an eigenvalue function 

Recall the dehnition of the sub differential to a convex function g : —)■ M at 77 ; 

dgiv) = {i^e I g{T) > g{ri) + : (r - 77 ) Vr G 

To receive the Mohr-Coulomb yield function or the plastic potential, we specify g as follows: 

^( 77 ) = an;i( 77 ) - bujsir]), 77 G (2.2) 
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where the parameters a,b > 0 are sufficiently chosen. Notice that the convexity of the eigenvalue 
function g can be derived from: 

a;i( 77 ) = maxr/: (e ® e) = max (r/e) ■ e, a; 3 ( 77 ) = min r/: (e ® e). (2.3) 

eeR3 eSK® 

|e|=l |e|=l |e|=l 


Specihc form of dg{r}) with respect to ( 2 . 2 ) can be found using a framework introduced in 


Chapter 2]. We derive another form of dg{r]) that is convenient for purposes of this paper. 
Lemma 2.1. Let g : —)■ M 6e defined by ( 2 . 2 ). Then for any r] G it holds: 

dgiv) = ji/ = ^ z/^ei ® e | (ei, ea, eg) eV{ri); a > z/g > z /2 > z /3 > - 6 ; 


2 = 1 


'^j^^ = a-b; - a)[ui{r]) - U2{r])] = 0 ; (z /3 + 6 )[(:u 2 (r/) - ^ 3 ( 77 )] = 0 I . (2.4) 


2 = 1 


Proof. Since 5 ^( 0 ) = 0 and g{2r]) = 2g{r]) the standard dehnition of dg{rf) is equivalent to: 

dgiv) = {i^ ^ I gip) = u :r]] ^(r) >u:t Vt G 
First, we derive necessary and sufficient conditions on z/ G ensuring 

g{T) >u:t Vt G 

To this end, consider the following spectral decomposition of u: 


V = 


^ vSi ® fj, Z/1 > Z/2 > Z^3, (fl, f2, fs) e V(ly) 


2 = 1 


Choose T — ±J, where I is the unit tensor in Then from (2.6), (2.7) we have: 

Vi + n.2 + 03 - a - b. 


(2.5) 

( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 


Choose T — fl ® f, and r — — fs ® fs. Then from (2.6), (2.7) we derive, respectively: 


Let T € 


d3x3 

^sym 


ui < a, z/g > —b. (2.9) 

be arbitrarily chosen and denote ti := t : (fj ® fj), i = 1, 2,3, (fi, f 2 , fs) G V{ly). Then, 

c<;i(t) > Ti > a; 3 (T), Vi = 1,2, 3, (2.10) 


Ti+T2 + T3 = T : I = UJi{t) + U2iT) + UJ3 {t), 
follow from I = Yl^=i f* ® f* (2-3), respectively. Consequently, 


ly : T = 


i^iTi = rfivi - z/2) + (ri + T2){v2 - + (ti + r2 + 


2=1 


- rfivi - z/ 2 ) + (r : / - Tfij{v2 - + lysT : I 

“ Wi(77)(z/i - z/2) + [r : / - ujfiri)]{iy2 - z/3) + z/gr : / = ^ 


2=1 


= i^i[c<;i( 77 ) - U 2 {'n)] + (z/i + z/ 2 )[a; 2 (r/) - ^ 3 ( 77 )] + (z/g + z /2 + iy3)oo3{'n) 


^ lyfiooiir]) -^2(77)] + {a-b- iy3)[(^2{'n) - i^ 3 (^)] + (a - b)uj3{r]) 




“ a[i:ui( 77 ) - ^ 2 ( 77 )] + a[i:u 2 ( 77 ) - ^ 3 ( 77 )] + (a - b)uj 3 {r]) 

= aufir) - bu3iT) = c/(r) Vr G R^y^. 


( 2 . 11 ) 
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Thus the conditions (2.7)-(2.9) are necessary and suffici ent for (2.6). 

Secondly, assume that v belongs to dg{rj). Then (2.7)-(2.9) hold. Since g{rj) u : rj, the 
equalities must hold within the derivation of (2.11) for t = rf, i.e., we have: 

(n - - z/2) = 0, (rg - a;3(r/))(z/2 - z/3) = 0, (2.12) 

- a)[u;i{r]) - uj 2 {'n)] = ^, {^^3 + b)^^) - (2-13) 


It is easy to see that the equalities in (2.12) imply: 


3 ( 61 , 62 , 63 ) € r ( 17 ) : 13 = ^l/j6.: 0 ej. 


(2.14) 


2 = 1 


We have proven that for any element G dg{r}) the conditions (2.7)-(2.9), (2.13) and (2.14) hold. 
Therefore, 


dg{ri) C < 


(g) e,; G 


d3x 3 
^sym 


( 61 , 62 , 63 ) G V{r]); a>Ui>U 2 >Us> -b; 


2=1 


= a - b; {ui - a)[ui{ri) - U 2 ir})] = 0; {^3 + b)[u 2 {ri) - usirf)] = 0 \ . (2.15) 


2 = 1 


Conversely, one can easily check that any element from the set on the right hand side in (2.15) 


belongs to dg{r}) using (2.5) and (2.11). 


□ 


Remark 2.1. One can easily specify the eigenvalues vi, z /2 and in (2.4) depending on a number 


of distinct eigenvalues of t}. If gi > ^2 > hs then z/i = a, z /2 = 0 and 1^3 = —b. If gi = r ]2 > ha then 
a > z/i > z/2 > 0, z/i + z/2 = a, and 1^3 = —b. If hi > h2 = hs then ui = a and 0 > z/2 > z/3 > —b, 
1^2 + 1^3 = -b. 


2.2 First and second derivatives of eigenvalue functions 

It is well-known that differentiability of eigenvalue functions depends on multiplicity of the eigenval¬ 


ues. For example, the function g is differentiable at rj with hi > h 2 > Vs as follows from Remark 2T 
Following [21 [22] , we derive the first and second Frechet derivatives of the eigenvalue functions using 

at 77 is denoted as VF{r]). Analogous 


3x3 

'"sym 




eigenprojections. The derivative of function F : 
notation, VF{r]), is also used for tensor-valued function F : —)■ Further, it is worth 

mentioning that some derivatives introduced below cannot be extended on 

First, assume three distinct eigenvalues of r/, i.e., hi > h 2 > hs- Then one can introduce the 
eigenprojections Ei := Ei{r]), i = 1,2,3, of r] as follows: 


F;* = 6i (g) 6i =-y— VlDip. —i = l,2,3. 
iVi - hi)(hi - hfc) 


It holds: 


V ^ ^ Vi^ii 'y 3 ^ 1 


2=1 


2=1 




V{rf 


Vui{r)) = Ei{vi), 7 = 1,2, 3, 

(Vj + hfc)I - (2hi - hi - Vk)Ei ®Ei- (vj - Vk) [Ej ® Ej - E^® E^] 


(hi - hi)(hi - Vk) 


(2.16) 

(2.17) 

(2.18) 
, (2.19) 
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for any i = 1 , 2 , 3 , i ^ j ^ k ^ i, where the components of the fonrth order tensors 'D{rf) and I satisfy 
\p{rf)]ijki = 5 ik['n\ij+ 5 ji[ri\ik and = 5 ik 5 ji, respectiveljQ We nse the notation Ei{r]) := VEi{r]), 
i = 1,2,3. 

Now, assnme f]! > 'r]2 > V3- this more general case, one can introdnce the derivatives of uj^ and 
uji2 '■= 1^1 +0J2- From (2.16), it is readily seen that the fnnction E^, can be continnonsly extended for 

( 2 . 20 ) 


T] satisfying t]i = r]2 nnlike Ei and E2- Hence and from (2.17), (2.18), one can write: 


Vuj^^ri) = Esir]), Vu^ir]) = I - E^{ri) =; Ei 2 {,'n). 

To continnonsly extend the fnnction ]E3(?7) := VE^ij]) = —VEi2{r]), we nse the eqnality 
(hi - h 2 )(^i ® El- E2® E2) = (r/ - hs^s) ® Ei 2 + Ei 2 ®{ri- i^sE^) - (r/i + 172)^^12 


E 


12 


and snbstitnte it into (2.19) for i = 3 . We obtain 

V{r]^) - {t]i + 772)1 - [r; ® E12 + E12 ® r/] + (771 + 772)^^12 


E3(r/) = 


E 


12 


(h 3 -hl)(h 3 -h 2 ) 

(771 + 772 — 2773)^^3 ® E^ + ? 73 [i^i 2 ® E2, + E^, ® Ei2\ 


+ 


(h 3 -hl)(h 3 -h 2 ) 


( 2 . 21 ) 


Clearly, (2.21) is well-defined also for 771 = 772. Notice that if 771 = 772 > 773 then rj has only 
two eigenprojections: ^^12 and E3, and 77 = r]iEi2 + rj^E^. Conversely, if 771 > 772 > 773, then 
Ei 2 = El E2. 

If 771 > 772 > 773 then one can introdnce the derivatives of the fnnctions Ui, UJ23 '■= 0J2 + i^3- 
Similarly as in the previons case, it holds: 


Vui{ri) = Ei{r]), Vu23{'n) = I - Eiir)) =: ^^23(^7), 


( 2 . 22 ) 


Ei(?7) = VEiij]) 


V{r]^) - (772 -F 773)1 - [?7 0 E23 + E23 ® 77] -F (772 -F r] 3 )E 23 ® ^23 

(hi -h2)(hi -h2) 


(h2 + h3 ~ 2771)^/1 ® El -|- rji [E23 ® El -f- El 0 E23] 
(hi -h2)(hl - h3) 


(2.23) 


Notice that if hi > h2 = hs then 77 has only two eigenprojections: Ei and E23, and 77 = 771^^1-)-773^^23- 
Conversely, if hi > h 2 > h35 then E23 = E2 + E3. 

In the general case hi ^ h2 > h35 h holds that hi + h2 + h3 = 'h • thns 


Vlui -t- a;2 -f 1^3] (77) = I. 


(2.24) 


Notice that if hi = h2 = h3 then 77 = rjil has only one eigenprojection: I. 

Remark 2.2. The mentioned derivatives can be fonnd in simpler forms when plane strain assnmp- 
tions are considered, see Appendix A of this paper. 


3 The Mohr-Coulomb constitutive problems 

In this section, we introdnce the Mohr-Conlomb constitntive initial valne problem and its implicit 
Enler discretization. We nse the model proposed in |2] containing the Mohr-Conlomb yield criterion, 
the nonassociative plastic flow rnle, and the nonlinear isotropic hardening. 

^In Appendix A], instead of 'D{rf') and I, their symmetric parts are introduced. For example, instead of I, the 
tensor Ig with the components [Is]ijki = + SuSjk) is considered. One can easily check that I : r; = Ig : r; = 77 

for any rj S A similar identity also holds for 'D{rf). 
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3.1 The initial value constitutive problem 


The initial value constitutive problem reads as: 

Given the history of the strain tensor e = e{t), t G [0,tmax]; the initial values £^(0) = £q, ^(0) = 
£o- Find (cr(t), ^(t)) such that 

eP€\dgM, = I (3.1) 

A > 0 , f{(T,K)<0, Xf{(T,K)=0. ) 
hold for each instant t G [ 0 ,tmax]- 

Here, cr, A denote the Cauchy stress tensor, the plastic strain, the hardening variable, and the 

plastic multiplier, respectively. The dot symbol means the pseudo-time derivative of a quantity. The 
functions / and g represent the yield function and the plastic potential for the Mohr-Coulomb model, 
respectively. They are dehned as follows: 


/(cr, k) = ( 1 -|-sin 0 )a;i(cr) — (1 — sin 0 )a; 3 (cr) — 2 (co-f k) cos 0 , 
g{cr) = ( 14 -sin'/)a;i(cr) - (1 - sin'/)a; 3 (cr). 


(3.2) 

(3.3) 


where ui and a ;3 are the maximal and minimal eigenvalue functions introduced in Section and 
he material parameters Cq > 0, 0, G ( 0 , 7 r/ 2 ) represent the initial cohesion, the friction angle, and 
the dilatancy angle, respectively. Notice that /, g are convex functions with respect to the stress 
variable. Recall that the function g was already introduced in Section 2.1 for the choice 


a := 1 -|-sin-/, b := 1 — sin'i/j 


(3.4) 


and thus one can dehne dgi^cr) using Lemma 2.1 

Further, the fourth order tensor Dg represents linear isotropic elastic law: 


Clearly, df{cr,K,)/dK = —2cos(/. 


cr = 


= ysA- 


2G){I : e^)I + 2Ge^ 


= 43A- 


2G)I ® / + 2GI, 


(3.5) 


where = e — is the elastic part of the strain tensor and K,G > 0 denotes the bulk, and shear 
moduli, respectively. 

Finally, we let the function H representing the non-linear isotropic hardening in an abstract 
form and assume that it is a nondecreasing, continuous, and piecewise smooth function satisfying 
H{0) = 0. 

It is worth mentioning that the value tmax need not be always known, see Section 


3.2 The discretized constitutive problem 


Let 0 = tQ < ti < ... < tk < ... < t]y = tmax be a partition of the interval [0,tmax] and denote 
CTfc := cr( 4 ), Sk := £(4), ^ := ^^(4), el ■= ^(4), := ^(4-i), 4'’ — ^(4) - £^(4-i), and 


err : = 


: . Here, the superscript tr is the standard notation for the so-called trial variables (see. 


e.g., [ 2 ]) which are known. If it is clear that the step k is hxed then we will omit the subscript k 


and write cr, s, and cr*^ to simplify the notation. The fc-th step of the incremental 

constitutive problem discretized by the implicit Euler method reads as: 

Given and Find cr, ^, and AA satisfying: 


cr = — AADe : jz, jz G dg{cr), 

= +AA(2cos4, 

AA > 0, f{a, H{^)) < 0, AXfia, H{A>)) = 0. 


(3.6) 


Unlike problem (3.1), the unknown is not introduced in (3.6). It can be simply computed from 


the formula s^{tk) = sifk) — Dg ^ : cr(ffc) and used as the input parameter for the next step. 
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4 Solution of the discretized constitutive problem 

The aim of this section is to derive an improved solntion scheme to problem 


The solntion 

scheme bnilds on the standard elastic predictor - plastic corrector method and its improvement is 


based on the form of dgi^cr) introdnced in Lemma 2.1 Within the elastic prediction, we assnme 
AA = 0. Then, it is readily seen that the triple 


cr = <T 


Ar 


eP = 


AA = 0 


is the solntion to (3.6) nnder the condition 

< 0 . 


(4.1) 


(4.2) 


The plastic correction happens when AA > 0. Then the nnknown generalized stress {cr,H{eP)) lies 
on the yield snrface and thus the corresponding plastic correction problem reads as: Given and 
A’’*''. Find cr, A’, and AA > 0 satisfying: 


cr = cr*^ — AADg : iz, jz G dg{cr), 
^ = +AA(2cos0), 


(4.3) 


/(<T,//(£!>)) =0. 


The solution scheme to problem (4.3) is usually called the implicit return-mapping scheme. Since its 


derivation is technically complicated, we divide the rest of this section into several subsections for 


easier orientation in the text. In Section 4.1, problem (4.3) is reduced and written in terms of principal 


stresses. In parallel Sections |4.2j|4.5| , we introduce solution schemes for returns to the smooth portion, 
to the “left” edge, to the “right” edge, and to the apex of the pyramidal yield surface, respectively. 


In Section 4.6, we derive a nonlinear equation for the unknown plastic multiplier. This equation is 
common for all types of the return and has the unique solution. Hence, we derive: existence and 


uniqueness of problems (3.6) and (4.3), a priori decision criteria for the return types, and other useful 


results describing a dependence of the unknown stress tensor on the trial stress tensor. 


4.1 Plastic correction problem in terms of principal stresses 

First, we reduce problem (4.3) using the spectral decomposition of cr (see Section]^: 


cr = (g) e*, Ui > (T2 > as, (ei, 62, 63) G I/(cr), ai := uji{cr), i = 1,2,3. 


2 = 1 


(4.4) 


From the dehnition of / introduced in Section it is easy to see that the equation ( 4 . 3)3 can 
be written only in terms the principal stresses cri,(T 2 ,cT 3 instead of the whole stress tensor cr. To 
re-formulate (4.3)i, we use Lemma 2.1 and ( 3.4[ ): there exists ( 01 , 62 , 63 ) G V(cr) such that u = 
Yfi=i Ae* ® Bi, where 


1 + sin-i/) > 1^1 > 1^2 > t'a > —1 + sin-i/), + 1^2 + = 2 sin' 0 , 

{ui -1 - sm'ijj){ai - (J2) = 0, (1/3 + 1 - sin'0)(cr2 - a^) = 0. 


(4.5) 


Since I = X)j=i ® (3-5) implies 


: jz = 


E 

2=1 


-{3K — 2G) sin-^ + 2Gvi 

3 


6i (g 6i. 


(4.6) 
























Then one can substitute (4.4) and (4.6) to (4.3 )i: 


cr*'' = cr + AADg : v = ® where = cr^ + AA 


2=1 


- {3K — 2G) sin tfj + 2Gui 
3 


. (4.7) 


Notice that (4.7)i dehnes the spectral decomposition of cr*'’. Since cxi > (72 > and > ^^2 > ^ 3 , 
we have: 

(*) at > at > at- 

{ii) if at = at then (j* = aj, Vi = Vj. 

From (i), it follows that the eigenvalues ati at, at are ordered and thus uniquely determined using 
the eigenvalue functions: at = utcr^^), i = 1,2, 3. From (ii), we conclude that cr = ® 

V = ® e*’’ for any (e^'’, e^'", 63 '’) G lA(cr*'’). The following lemma summarizes the proven 

results. 


Lemma 4.1. Let {(j,^,/W), AA > 0 6 e a solution to (4-3) for given and A’’*'’. Let ai, at, 
i = 1 , 2 , 3, he the ordered eigenvalues of cr and cr*^, respectively. Then {ai,a 2 , a^, A', AA) is a solution 
to: 


ai = a, 


tr 


A A [|(3iF — 2G) sin'^ + 2Gof\ , i = 1, 2, 3, 
^ = +AA(2cos0), 

(1 + sin0)(7i — (1 — sin(/))(73 — 2(co + H{^)) coscj) = 0, 


(4.8) 


where i^i,r' 2 A 3 satisfy (f-S) Conversely, if {ai,a 2 ,a^,e^, A\), AA > 0 is a solution to (4-8) then 
(cr,A’, AA) solves (4-3), where cr = ® ^ 2 '’, 63 ’’) G lA(cr*'’). 


To be in accordance with problems (3.6) and (4.3), we do not include h'i,h' 2 ,r '3 to the list of 


unknowns. From (4.5), it follows that the values of z/i, z/ 2 , 1^3 can be specihed depending on multiplicity 


of ( 7 i, ( 72 , ( 73 , similarly as in Remark |2.1[ Therefore, we will distinguish below four types of the return 
on the yield surface: the return to the smooth portion fyi > a 2 > ( 73 ), the return to the left 
edge (( 7 i = a 2 > ( 73 ), the return to the right edge {ai > a 2 = ( 73 ) and the return to the apex 
(( 7 i = (72 = ( 73 ). This terminology follows from |2], another one is used, e.g., in [5]. Within the 
below introduced notation, we will use the subscripts s, I, r, a to distinguish the return type and the 
superscript “tr” to emphasize a known quantity depending only on the trial variables. 


4.2 The return to the smooth portion 


Assume ai > a 2 > ( 73 . Then vi = 1 + sin-^, z /2 = 0 , z /3 = —(1 — sin'^) and (4. 8)1 reads as: 


ai = at-AX 


(72 = at — AA 


^3 = at-AX 


-{3K — 2G) sin-^ + 2G{1 + sin-^) 
3 


-{3K — 2G) sin-^ 
3 


-{3K — 2G) sin-^ — 2G(1 — sin if) 
3 


(4.9) 

(4.10) 

(4.11) 


Consequently, one can substitute (4.9), (4.11), and ( 4 . 8)2 to ( 4 . 8 ) 3 . This leads to the equation 
gf (AA) = 0 , where 


(iTil) = (1 + sin(/))(7f — (1 — sin(/))(73'’ — 2 [co + FT (e**’*” + 7(2cos 0))] cos(/) 


-7 


- {3K — 2G) sin if sin cf + 4G(1 + sin if sin 0) 
3 


(4.12) 


Further, from (4.9)-(4.11), two additional important consequence follow: 
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(4.13) 


af > a*" > a*", 


AA e Af := {7 e (0, +00) | 7 < min{7*0,7‘(J}, where 


rr^r _ tr 

■■= .Z\ , > 0 > 71 :.: = 


2G{1 + sin-i/)) 


rr^r _ tr 


> 0 . 


4.3 The return to the left edge 

Assume ai = a 2 > cs. Then 1/3 = —(1 — sin-^), ui + U 2 = 1 + sin-i/), and 1 + sin-^ > i^i > U 2 > 0 


implying ui — U 2 < I + sin-i/). Consequently, (|4.8|)i yields: 

1 


+ = cri = + 


as = aJ — AA 


AA 


-{SK — 2G) sin'^ + (7(1 + sin'^) 

o 


-{3K — 2G) sin'^ — 2(7(1 — sin-^) 
3 


and 


0 = cTi — ct 2 = (jf — a 2 ~ AA[2(7(i/i — i> 2 )] > erf — erf — AA[2(7(1 + siu'^) 
After substitution (4.14), (4.15), and ( 4 . 8)2 to ( 4 . 8 ) 3 , we arrive at gf (AA) = 0, where 


(4.14) 

(4.15) 

(4.16) 


sin 0 )(crf + erf) - (1 - sin 0 )crf - 2 [cq + i7 + 7 ( 2 cos 0 ))] cos 0 


7 


-{3K — 2(7) sin'0sin 0 + (7(1 + sin'0)(l + sin0) + 2(7(1 — sin'0)(l — sin0) 
3 


(4.17) 


Further, from (4.14)-(4.17), three additional important consequences follow: 


^f > 4 ", 


eri, as, A A depend on erf, erf only through erf + erf, 
AA e Gf := {7 e (0, +C50) | 7 ^ < 7 < 7 f„}, where 
erf + erf — 2erf 1 + sin 


fv 

%a = 


tr , / 1 l + simjj . 


2(7(3 —sin fi) 3 —sin'^ 


7fz + 1 


3 — sin ^|J 


7„ > 0 


(4.18) 


and 7 f, 7 :f are the same as in (4.13). Notice that 7 f < < 7 f^ in this case. 


4.4 The return to the right edge 

Assume eri > er 2 = as- Then Vi = 1 + sin-^, 1^2 + Vs = + sin-e/i, and 0 > z /2 > ^^3 > —1 + sin-e/i 


implying Z 4’2 — z /3 < 1 — sin-^. Consequently, (4.8 )i yields: 


eri = erf-AA 
1 


-{3K — 2(7) sin'^ + 2(7(1 + sin'^) 
3 


-(er2 + er3) — er3 — ^(erf + erf^ 


AA 


-(3A' — 2(7) sinfi — (7(1 — sin'^) 

o 


(4.19) 

(4.20) 


and 


0 = er2 — er3 = erf — erf — AA[2(7(z/2 — 1^3)] > erf — erf — AA[2(7(1 — sin^’)]- (4.21) 
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After substitution (4.19), (4.20), and ( 4 . 8)2 into ( 4 . 8 ) 3 , we arrive at g*''(AA) = 0, where 


(iTil) = (1 + sin(/))crf - -(1 - sin 0 )(cr 2 '’ + Ug'’) - 2 [co + H + 7 ( 2 cos 0 ))] cos 0 - 


7 


-{3K — 2G) sinsin 0 + 2G(1 + sin0)(l + sin0) + G(1 — sin-0)(l — sin0) 
o 


. (4.22) 


Further, from (4.19)-(4.22), three additional important consequences follow: 

. af > ar > 

• (Ti, ( 73 , AA depend on af only through , 

• AA e := {7 e ( 0 , +oo) | < 7 < 7 *(„}, where 

2a^{ — 1 — sin 0 


7*’’ = 

ir.a 


2G{3 + sin 0) 3 + sin 0 


7; + (1 _ FiftT; 7-> 0, 

1 3 Asin^ ' “ 


(4.23) 


and 7 *'), 7 *)), are the same as in (4.13). Notice that 7 *)), < 7 *^ < 7 *') in this case. 


4.5 The return to the apex 

Assume ai = a2 = (73. Then ui + i>2 + I's = 2 sin '0 and 1 + sin^ > > 1^2 > ^3 > — 1 + sin^ 


implying 2 z/i — z/2 — z/3 < 3 + sin^, ui + U2 — 2 v^ < 3 — sin^. Consequently, ( 4.8 )i yields: 
(7i = -{ci + a2 + cr^) = + ^^2 + ^3^) ~ AA[2iF sin-i/’] 

and 


0 = 2(7i — (72 — (73 > 2(7^ 


tr 


Ar 


(To 


Ar 


(To 


AA[2G(3 + sin: 


0 = (7i + (72 — 2(73 > + (^2 ~ 2(73 ~ AA[2G(3 — sin0)] 


(4.24) 


(4.25) 

(4.26) 


After substitution (4.24) and ( 4 . 8)2 into ( 4 . 8 ) 3 , we arrive at g*''(AA) = 0, where 


( 7 ) = + ^2 + sin0 — 2 [co + if (A’’*'" + 7(2 cos0))] cos0 — 7 [ 4 iF sin0 sin 0]. (4.27) 

o 


Further, from (4.24)-(4.27), two additional important consequences follow: 

• ( 7 i, AA depend on af, a^i only through af + 

• AA e C'*’’, := {7 G (0, +CX)) I 7 > max{ 7 ;* 0 , 7 * 0 }}, where 7 }*}, Ar,a same as in 

(|4.18l), (|4.23l), respectively. 


4.6 Solvability analysis and a priori decision criteria 


In parallel Sections |4.2j|4.5| , the solution schemes for the investigated return types were introduced. 
Similar schemes are also known from literature (see, e.g., [21 Section 8 ]) and their solutions are 


candidates on the solution to problem (4.8). This current approach is based on a blind guesswork 


since the position of the stress tensor on the yield surface is not a priori known. However at the 
ends of Sections [4.2 4.5, we also derived some additional results following from (4.5), i.e., from the 
knowledge of dg{(T). These results enable to improve the solution scheme to problem (4.8). First, 
we use the sets (7*^, Gf^, C}'’, G^, the values 7 *)), 7 * 0 , Ai% 7^0; die equations qf{AX) = 0, 
qj^{AX) = 0, g0(AA) = 0, g0(AA) = 0 introduced above to find a unique nonlinear equation for the 
unknown plastic multiplier. 
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Lemma 4.2. There exists a unique function : M_|_ —)■ M satisfying: 


/■\ ^tr\ _ ^tr ^tr\ _ ^tr r^trl _ ^tr ^tr\ _ ^tr 

[l) Q Icf - Qs ’ <1 \cT - (ll ’ Q lc‘- - Qr ’ <1 \ci- - <la ■ 

{ii) is continuous, piecewise smooth and decreasing in 

(m) g*^(0) = /(cr*^i7(^’*^)). 

{iv) q^^i'y) —)■ —oo as 'y ^ +cxd. 


Proof. Notice that the values 7 *^;, 7 *^^, and 7 *^^ are nonnegative and a priori known. Moreover, 
from (4.13), (4.18) and (4.23), it follows that only two ordering of these values are possible: either 

iZ < iZa < <a < itr Cr < < Zip 

First, assume < 7*^^. Then Cl = (0,7^yj, = [Zip if, I^ ^'1 = Ka^+oo)- 

Dehne the function 


g*^(7) = (l + sin0)crr-(l-sin0)crr-7 


-{3K — 2G) sinsin0 + 4G(1 + sin0 sin0) 
3 


+G(1 + sin0)(l + sin 0)(7 - 7 ,C 0 + + -G(3 - sin0)(3 - sin 0)(7 - 7 ^;)+ 

-2 [co + + 7(2cos0))] COS0, 76 ( 0 ,+cx)), (4.28) 

where (.)■*■ denotes a positive part of a function. It is easy to verify that g*^ has the required properties 
under the assumptions on H from Section 

Secondly, assume 7 *^^ < 7 *y;. Then (Ff = ( 0 , 7 *yj, = 0, = ZlZZl,!^ = [7r(a)+oo) and 

the function g*'’ with the required properties is dehned as: 


g^( 7 ) := ( 1 + sin0X - ( 1 -sin 0 )a 3^-7 


-{3K — 2G) sin0 sin 0 + AG{1 + sin0 sin0) 
3 


+G(1 - sinV^)(l - sin0)(7 - + -G(3 + sinV^)(3 + sm(j)){-f - 7 ^yj 

— 2 [co + H + 7(2 COS0))] COS0. 


(4.29) 

□ 


Remark 4.1. Notice that formulas (4.28) and (4.29) coincide for 'jZ = 'yll- Hence, the function 
g( 7 ; (T 0 , CT 0 , (T 0 , = Z^iZ is continuous and piecewise smooth with respect to the trial variables. 


Lemmas 4.2, 4.1 and (4.2), (4.1) imply the following main results. 


Theorem 4.1. Let g*'’( 0 ) = f{cr^^,H{e^'Z) > 0 . Then the equation g*''(AA) = 0 has a unique 
solution in M+. The solution vanishes if and only if f , H{e^'Z) = 0 . Moreover, if there are 
71,72 > 0 such that 71 < 72 , g*'’( 7 i) > 0 , and g*^(72) < 0 , then AA G (71,72)- 


Theorem 4.2. Let f{(T^'^,H{e^’Z) > 0. Then problems (4-3) and (f.S) have a unique 
solution components to problem (4-8) can be found in the following way: 


solution. 1 fie 


1. Let g0(min{7*0 7*0}) < 0. Then AA G Cf is the unique solution to g0(AA) = 0 and ui > 


2 . 


<72 > <73 can be computed from (4.9)-(4.11). Moreover, a} 


> 


> a. 


Let Zh 


< 


ZZ Znziz > 0 and qZiZZ < 0. Then AA G 


cr 


g0(AA) = 0 and cxi = ct 2 > <73 can be computed from (4-14h and (4-15). Moreover, > a\ 


and AX, ai, as depend on aZ, aZ only through cr0 + aZ 


is the unique solution to 

tr 
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3. Let Qti.'lXr) ^ 0 < 0- Then AA G is i/ie unique solution to 

g*^(AA) = 0 and cxi > ct 2 = cxs can 6e computed from (4-19), and (4-20). Moreover, a\^ > 
and AA, cxi, as depend on a^, a^ only through a^ + al^. 

4- Let g*’'(max{ 7 *’^, 7 *^^}) > 0. Then AA G is the unique solution to q*J{AX) = 0 and 
(Ti = ct 2 = as can be computed from (4-24)■ Moreover, AA and ai depend on af, a^ and a^ 
only through aff + a^{ + ■ 


The component A’ can be computed from (4-8) for all of these cases. 


Theorem 4.3. The discretized constitutive problem (3.6) has a unique solution. 


Remark 4.2. Notice that Theorem 4.2 contains the solntion scheme to problem (4.8) and snmmarizes 
the advantages of the snbdifferential treatment within the constitntive solution: 

1. Existence and uniqueness of the solution. This expected result is not usually discussed in 
literature. 


2. A priori known decision criteria. Such criteria were known only for linear function H (see, 
e.g., [5]) where the solution components can be found in closed forms. 

3. Dependence of ai,a 2 , as, AA on af", alf, a^4f has been described in more detail than it is known 
from literature. This enables us to simplify construction of the stress-strain and consistent 
tangent operators introduced in the next section, and discuss semismoothness of the stress- 
strain operator. 


5 Stress-strain and consistent tangent operators 


In this section, we extend the solution scheme from Theorem |4.2| to problem ( |3.6[ ) and dehne the 
stress-strain operator and its derivative, i.e., the consistent tangent operator. Beside the results from 
Section]^ we also use the framework from Section [2^ based on eigenprojections and their derivatives. 
The stress-strain relation can be represented by an implicit function T; 


cr(4) ;= T {e{tk)-, £p(4-i), A’(4-i)) • 


If we £x step k and recall = D^itk-i), = s(tk) — s^{tk-i), one can write 


o- := T(s;£P(4-i),£^(4-i)) = 5 


(5.1) 


omitting the subscript k. The consistent tangent operator for step k will be represented the Frechet 
derivative VS = V^trS. If it exists at (£*'’, A’’*'’) then V^T = V^trS. It is sufficient to derive the 
operators S and VS on the following open sets: 


Mf 

Mf 

Mr 

Mr 

M^^ 


Ar 


^tr 


^tr 


Ar 


^tr 


d3x3 

^sym 

d3x3 

^sym 

[)3x3 

^sym 

[)3x3 

^sym 

d3x3 

^sym 


J-O \ / 


gf(0) > 0, gf(min{7^(;,7^yj) < 0}, 

7S < Tta. > 0, qtr-fta) < 0}, 
< < <a, rrht) > 0, gr«a) < o}, 

tr / r i.r f.r x \ ^ 


g^'’(max 


HaMra}) > 0} 


From Section]^ it follows that these sets are mutually disjoint and the closure of their union is equal 
to since cr*'’ = Dg : Further, the tensors cr*'’ and have the same eigenvectors and their 

eigenvalues are related as follows: 


a. 


tr 


= -{3K- 2G){er + 4 -1- 4") + 2Ge 
3 


i = 1,2,3. 


(5.2) 
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Hence, ef > if and only if for any i,j = 1, 2, 3. Therefore, cr*'" and also have the same 

eigenprojections. For the sake of simplicity, we assume that H is differentiable at + AA(2cos0) 
and denote Hi := + AA(2cos0)). 

The elastic response. Let G . Then, clearly, 

S = De : VS (e'", = De. (5.3) 


The return to the smooth portion. Let e*'’ G M*'’. Then > 62 > ^ 3 ^ holds and consequently, 
the values := Ef := Efe^’’), i = 1,2,3, are well-dehned as follows from Section 2.2 

Therefore, 

3 3 

(5.4) 



5(e 

tr 

= Y.a,Ef, P5(s*fe^’’^ 

T = E 

[(jiE 

f + Ef ® Vaf\ 


Since, 










Vai 


Im- 

2G)I + 2GEf - 

■T>(AA) 

Im- 

2G) 

sin'0 

+ 2G(1 + 

sin'0) 

Va2 

pToli 

Im- 

2G)I + 2GEf - 

■T>(AA) 

fiK- 

2G) 

sin'0 



Va^ 

ED 

Im- 

2G)I + 2GEf - 

■T>(AA) 

fiK- 

2G) 

sin'0 

-2G(1- 

sin-^) 

we have 










VS 

■) = 

3 

f + 2GEf ® Ef 

] +^(3iL-2G)i 

®I 

— 




2=1 


2G(1 + sin - 2G(1 - sin V’)-Ef + -(3iF - 2G) sinV^J 


T>(AA), (5.5) 


where 


/A 1 S 2G(l + sin0)F;f-2G(l-sin(^)F;f+ |(3iL-2G')sin(^/ 
I (3iL — 2G) sin xjj sin 0 + 4G(1 + sin i/j sin 0) + 4iLi cos^ 0 


The return to the left edge. Let G Mf. Then ai = (J 2 and ef > ef. It means that one can 
introduce the notation := Ef ’2 .= I — Ef, Ef := E 3 (e*'’) and write 

S (£*f = aiF;f2 + ^3-Ef, VS (£*f e^’*") = (ag - cTi)Ef + Ef ® Vai + Ef 0 Va^. (5.6) 

Since, 


Vai ^ hsK-2G)I+ GEf2-V{AX) 

o 

Vas ® ^(3iL-2G)/ + 2G'F;f-T>(AA) 


- (3iL — 2G) sin V’ + 1^(1 + sin fj) 
o 


-{3K — 2G) sinf> — 2G(1 — sin-^) 
3 


we have 


VS (e'f eP’'") = (ag - afEf + GEf 2 ® Ef 2 + ‘iGEf ® Ef + - {3K - 2G)I ® / - 


G(1 + sinf>)F;f 2 - 2G(1 - sini))Ef +-{3K - 2G) sinf>/ 


P(AA), (5.7) 
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where 


V{AX) ^ 


^(1 + sm^)E% - 2G{1 - sin (j))El^ + 1{3K - 2G) sin0/ 

I {3K — 2G) sin -0 sin 0 + G(1 + sin -0) (1 + sin (j)) + 2G(1 — sin -0) (1 — sin 0) + AHi cos^ (f) 


The return to the right edge. Let G M*'’. Then ct 2 = cxs and ef > e ^2 ■ ii' means that one can 
introduce the notation E^{ ;= Ei{e^^)., E 23 ^ ~ ^*1 1 •= IEi(£*’') and write 

S (£'^ = aiE^{ + a^E%, VS = (ai - a3)Ef + E^{ ® Vai + E% ® Va^. ( 5 . 8 ) 

Since, 


Vai ® J(3iL-2G)/ + 2G^;f-T>(AA) 
Va3 ® J(3iL-2G)/ + G^;*3-P(AA) 

o 


-{3K — 2G) sin'0 + 2G(1 + sin'0) 
3 


-{3K — 2G) sin-^ — G(1 — sin'0) 
3 


we have 


VS A’’'") = (cTi - a3)Ef + 2GE{^ ® E\^ + GE^^^ ® E% + -{3K-2G)I®I 

2G(1 + sinV')^;f - G(1 - sin?/^)^;^)^ + -{3K- 2G) sini/'J 


P(AA), (5.9) 


where 


T>(AA) “ 


2G(1 + sin (t))E^{ - G(1 - sin (i))E% + f (3iL - 2G) sin^J 
I {3K — 2G) sin -0 sin cj) + 2G(1 + sin -0) (1 + sin 0) + G(1 — sin ■0) (1 — sin 0) + 4ifi cos^ (f) 


The return to the apex. Let G M*'”. Then ai = a 2 = cr^ := p and 


S =pl, p = p'" - {2Ksinij)AX, p'" = -« + a*'' + a*") = K{e\^ + A/ + 4"), (5.10) 

vs (e-. eO-'-) S a AV 0 / = (1 - 2 /f sm /f / » /. 

Here, we use = KI. From the implicit equation g*^(AA) = 0, we obtain 


dAX JTWl 


sm( 


dp^^ 2K sin A sin 0 + 2Hi cos^ 0 


Hence, 


VS 


=k(^i 


K sin A sin 0 


K sin A sin (j) + Hi cos^ 0 


I® I. 


( 6 . 11 ) 


Remark 5.1. For each of the return type, we derived just one formula for VS without any other 
branching that depends on multiplicity of e^{,e^ 2 A^ 3 - This was achieved due to deeper analysis of 
dependencies within the constitutive solution, see Theorem 4.2 The additional branching in VS is 


introduced, e.g., in [21 Appendix A]. In many other references, VS is correctly derived only under the 
assumption > 8^2 > ^ 3 - However, such formulas can cause signihcant rounding errors in vicinity 
of the multiple eigenvalues. 
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Remark 5.2 

, 

of T 




Notice that one can continuously extend the dehnition of T (■; £^(ffc_i), = 

U M*'' U M^). Further, one can investigate semismoothness 


R™ \ (Vf U Mf U Ml 


on 

;£P(tfc_i),£^(tfc_i)) in This property ensures superlinear convergence of algorithms 

introduced in the next section. To show the semismoothness in M*'', M/'’, and one 

can use a standard framework introduced, e.g. in dHl EH 123 [HI EDI]. At the remaining points, the 
semismoothness is also expected based on Remarks |4.1 and 4.2 j 


However, its eventual proof seems 


to be more involved and we will skip it for the sake of brevity. 


Below, we use the notation T (■; {tk-i), ^{tk-i)) for the Clark generalized derivative of T with 
respect to the strain tensor. Clearly, T = VsT {s;e^{tk-i),^itk-i)) when 

T (•; £^(tfc_i), is differentiable at e. 


6 Direct and indirect methods of incremental limit analysis 

Inserting the stress-strain operator T to the balance equation, we obtain the incremental boundary 
value elastoplastic problem lan. This problem is further discretized in space by the hnite element 
method and combined with the limit load analysis as it is usual in Mohr-Coulomb plasticity 121 ED . 
In this section, we introduce the direct and indirect methods of the incremental limit analysis. For 
the sake of brevity, we focus only on an algebraic formulation of the problem. 

The vector of internal forces and the consistent tangent stiffness matrix at the fc-th step are 
represented by functions Fk : M"" —)• and Kk : M” —)■ respectively. It is worth mentioning 

that Fk and Kk are assembled using the operators T and T at each integration point [D- Notice that 
the algebraic representation of the used second and fourth order tensors is introduced in Appendix 
B. Further, we consider the load of external forces at step k in the form C,kl where Z G M"" is hxed 
and Cfc := C(^fc). Then the fc-step problem reads as: 

(T’fc)c given C,k e M+, hnd UkeW^ ■. F k{uk) = (kl, 

where Uk is the displacement vector. We assume that the parameters C, and t coincide and their 
limit value (um is unknown. It is well known that the investigated body collapses when this critical 
(limit) value is exceeded. Therefore, (um is an important safety parameter and beyond Czim no 
solution exists. Possibly Czim = +oo, however in meaningful settings of the problem, (nm is hnite. 
The simplest computational technique is based on the so-called incremental limit analysis where we 
adaptively construct the sequence 


0 < Cl < C 2 < • • ■ < Cfc < Cfc+l <...<C,lim 

depending on solvability of {Vk)c_ to detect inadmissible load factors. In practice, the increment of Cfc 
decreases when a chosen numerical method does not converge at step k. Such blind determination 
of Czs is an evident drawback of this direct method. 

More sophisticated adaptive strategy is based on local and/or global material response of the body 
on the prescribed load history. To this end, we compute the values ak = b^Uk, k = 1,2,... where Uk 
is the solution to (Vk)( and b is chosen so that to be the sequence { 0 ^} increasing. There are many 
ways how to do it. For example, one can detect a point on the investigated body where it is expected 
that a selected displacement is the most sensitive on the applied forces. Then b is the restriction 
of the displacement vector to its component. More universally, one can also set b = 1. This choice 
represents the work of external forces, is meaningful even for continuous setting of the problem and 
was analyzed in [251 EH ED EB] for generalized Hencky’s plasticity. Clearly, if the increment ak — ak-i 
signihcantly enlarges with increasing k then it is convenient to reduce the increment of C, for the next 
step. 
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The knowledge of suitable b also enables to introduce the indirect method of incremental limit 
analysis where the increasing sequence [ak] is given and the sequences {Ca:} and {uk} are computed 
using the following auxilliary problem: 


(■Pfc)" given (h, ak) G M" x M+, find {uk, (k) G x M+ : 


b'^Uk = ak- 


Clearly, if (ukXk) is the solution to ("Pfc)" then Uk also solves (Vk)( for X and (k < Cam < +C) 0 . 
Unlike to problem {Vk)c,^ one can expect that problem has the solution for any ak- Since the 

parameter a can be enlarged arbitrary, the indirect method is more stable and does not include any 
blind guesswork unlike the direct one. This is the main advantage of the indirect method. For the 
associative Mohr-Coulomb model, one can expect that C,k —t Cum as ak —)■ +cxd. This is proven in 
PUII7F] ioi b = I and the generalized Hencky’s plasticity. For the nonassociative Mohr-Coulomb 
model with ip « p, we observe that X ~ Cum for some hnite k and for k > k, the sequence {Cfc} is 
nonincreasing. In such a case, the material exhibits softening behavior and the direct method is too 
convenient. It is also worth mentioning that the indirect method is similar to the arc-length method 
introduced, e.g., in [23112]. 

We solve problems {V)t^ and ("P)" by the semismooth Newton method: 

Algorithm 1 (ALG-C). 

1: initialization: 

2: for z = 0,1, 2,... do 

3: find^w^GV: Kk{'u}^)5X = C,kl - Fk{K) 

4: compute = u\ + 6u^ 

5: if ||(5it*||/(||u*+^|| llit^ll) < CNewton then stop 

6: end for 
7 : set Uk = 


Algorithm 2 (ALG-a). 

1: initialization: 

2: for z = 0,1, 2,... do 

3: find v\ re* G V: Kk{u\)v'^ = Qd - Fk{u\), Kk{u\)w^ = I 

4 : compute 5Q = [ak — b^ 

5 : compute + SCpw'' 

6: set = K + du\ = Cfc + K 

7: if ||(5n*||/(||n^+^|| -F ||n*^||) < CNewton then stop 

8: end for 

9: set Uk = C,k = C^^- 


If T (•; £P(ffc_i), P’(ffc_i)) is semismoothness in then one can easily show that Fk is semis¬ 
mooth in M". The semismoothness is an essential assumption ensuring local superlinear convergence 
of these algorithms (see, e.g., [26]). Further, we initialize ALG-C and ALG-a using the linear extrap¬ 
olation of the solutions from two previous steps. In particular, we prescribe 

Uf. = Uk-i H- [Uk-i - Uk-2), Ck = Ck-I H-(Cfc-i - Ck- 2 ) 

«fc-l — (y-k-2 ctfc-l — Ctk-2 
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in ALG-a for k > 2, and analogously, in ALG-C- We observe that this initialization is more convenient 
than ul = Uk-i, Cfc = Cfc-i- 


The direct and indirect methods of incremental limit analysis are compared in Section 7.1 


7 Numerical experiments - slope stability 

We have implemented the direct and indirect methods of incremental limit analysis in MatLab for 
3D slope stability problem and its plane strain reduction. These experimental codes denoted as 
SS-MG-NP-3D, SS-MG-NH and SS-MG-NH-Acontrol are available in [30]. The codes are vectorized 
and include the improved return-mapping scheme for the Mohr-Goulomb model in combination with 
ALG-C or ALG-a. One can choose: a) several types of finite elements with appropriate numerical 
quadratures; b) locally rehned meshes with various densities. 



Figure 1: Cross section of the body with the coarsest mesh for Q2 elements. 

We consider the benchmark plane strain problem introduced in [2], Page 351] and its extension 
for 3D case. The 2D cross-section of the body with the coarsest mesh considered in [301 SS-MG- 
NH] is depicted in Figure The 3D geometry and the corresponding hexahedral mesh arise from 
2D by extruding. The slope height is 10 m and its inclination is 45°. On the bottom, we assume 
that the body is fixed and, on the lateral sides, zero normal displacements are prescribed. The 
body is subjected to self-weight. We set the specific weight pg = 20kN/m^ with p being the mass 
density and g the gravitational acceleration. Such a volume force is multiplied by the load factor 
C. The parameter a is here the settlement at the corner point A on the top of the slope to be in 
accordance with [2|. Further, we set E = 20 000kPa, u = 0.49, (j) = 20° and c = 50kPa, where c 
denotes the cohesion for the perfect plastic model. Hence, G = 67114kPa and K = 3 333 333 kPa. 
The remaining parameters of the Mohr-Goulomb model will be introduced below depending on a 
particular experiment. 

We introduce one experiment for the plane strain (2D) problem and one for the 3D problem. The 
primary aim of these experiments is to numerically illustrate that the formulas derived in Sections 
and Appendix A work well. This can be confirmed by observing the superlinear convergence of 
ALG-C and ALG-a and their stability in vicinity of the limit load. We also prescribe a high precision 
of these algorithms by the setting ej^ewton = 10“^^ in both experiments. Other aims will be specified 
below. 

7.1 Comparison of the direct and indirect methods in 2D 

We compare the direct method (code SS-MG-NH) and the indirect method (code SS-MG-NH- 
Acontrol) of the incremental limit analysis on the slope stability benchmark in 2D. We consider 
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the associative Mohr-Coulomb model containing the nonlinear isotropic hardening dehned as in pQ: 




min 


C — Co, 


- 


4(c - Co) 



Co = 40 kPa, H = 10000 kPa. 


Here, H represents the initial slope of H and the material response is perfect plastic for sufficiently 
large values of The function H is smooth and its influence on the limit load factor is negligible 
based on expertise introduced in [1]. We set V’ = 0 to have the associative model. 

Further, we use the Q2 elements (i.e. eight-noded quadrilaterals) with 3x3 integration quadrature 
and the mesh with 37265 nodal points including the midpoints and with 110592 integration points. 
The mesh has a similar scheme as in Figure [^but, of course, it is much more hner. Since the Matlab 
code is vectorized, we £x 10 inner Newton’s iterations for finding the unknown plastic multipliers in 
each integration point. 

Recall that in each step of the direct method, we solve problem {V)c_ using ALG-C. We set 
the initial load increment 5C,o = 0.5. If ALG-C converges during 50 iterations for step k > 1 and 
if the computed increment of the settlement satishes a*, — ak-i < 0.5 m then we set ^Cfc+i = ^Ck- 
Otherwise, the increment is divided by two. Within the indirect method where problem (P)" is 
solved using ALG-a we set the initial increment 6ao = 0.0414 of the settlement to have comparable 
results with the direct method. If the computed load increment satishes \(k — Cfc-il > 5e — 3 then we 
set 6ak+i = Sak- Otherwise, Sak+i = 25ak- The loading process is terminated when the computed 
settlement exceeds 4 meters for both methods. 




Figure 2: Load path for the direct method. Figure 3: Load path for the indirect method. 




Figure 4: Number of iterations for ALG-^. Figure 5: Number of iterations for ALG-a. 
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Figure 6: Convergence at selected steps for the Figure 7: Convergence at selected steps for the 
direct method. indirect method. 


The comparison of the direct and indirect methods is depicted in Figures [2]j^ The resulting 
loading paths practically coincide for both methods and they are in accordance with (am EH]. The 
computed limit value is equal to 4.057 which is close to the estimate 4.045 known from |23]. Both 
methods need 18 load step and have superlinear convergence in each step. Their convergence is 
similar up to step 11. However, other comparisons turn out that the indirect method behaves better 
than the direct one. First, the indirect method has less number of iterations between steps 12 and 
18. Secondly, the direct method contained 8 additional load steps without successful convergence 
while the indirect one convergences in each step. The successful load steps for both methods are 
depicted by the circular points in Figures and respectively. We see that the positions of these 
points are more convenient in Figure]^ than in Figurewith respect to the curvature of the loading 
path. Thirdly, we see in Figure]^ that the convergence in steps 11 and 16 is superlinear only up to 
le-10. Then, values of the stopping criterion oscillate. This is also observed for a few other steps 
of the direct method (e.g., steps 14 and 15). For the indirect method, this is not observed at any 
step. Finally, the computational times of the direct and indirect methods on a current laptop were 
approximately 9 and 7 minutes, respectively. 


7.2 Associative perfect plastic 3D problem 


Within the 3D slope stability experiment (code SS-MC-NP-3D), we compare the loading paths for 
the Q1 and Q2 hexahedral elements with 8 and 20 nodes, respectively. We consider 2x2x2 and 
3x3x3 noded integration quadratures for these element types, respectively. Two hexahedral meshes 
are prepared for this experiment. For the Q1 elements, the meshes contain 5103 and 37597 nodal 
points, 34560 and 276480 integration points, respectively. For the Q2 elements, the meshes contain 
19581 and 147257 nodal points, 116640 and 933120 integration points, respectively. We use the direct 
method of the incremental limit analysis which is terminated when the computed settlement exceeds 
5 meters. 

The corresponding loading paths are depicted in Figure We observe that the estimated limit 
values of ( are close to the expected value of 4.045 for the Q2 elements but not for the Q1 elements. 
To estimate (um using the Q1 elements, it would be necessary to use much hner meshes. Figures]^ 


and 10 illustrate failure at the end of the loading process for the Q2 elements and the hner mesh. 


8 Conclusion 

This paper extended the subdifferential-based constitutive solution technique proposed in [T] to 
elastoplastic models containing the Mohr-Coulomb yield criterion. It enabled deeper analysis of 
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Figure 8: Comparison of the loading paths for Q1 and Q2 elements. 



Figure 9: Total displacement and deformed shape Figure 10: Plastic multipliers at the end of the 
at the end of the loading process. loading process. 
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the constitutive problem discretized by the implicit Euler method and consequently led to several 
improvements within solution schemes. For example, a priori decision criteria characterizing each 
type of the return-mapping were derived even if the solution could not be found in closed form. The 
construction of the consistent tangent operator was also simplihed. Moreover, the paper brought self- 
contained derivation of the constitutive operators which is not too often in Mohr-Coulomb plasticity. 

The improved constitutive solution schemes were implemented within slope stability problems in 
2D and 3D. To this end, the direct and also the indirect methods of the incremental limit analysis 
were used and combined with the semismooth Newton method. Local superlinear convergence in 
each step of both methods was observed. Further, it was illustrated that the indirect method led 
to more stable control of the loading process or that higher order hnite elements reduced strong 
dependence on mesh. 


Acknowledgements 

The authors would like to thank to Pavel Marsalek for generating the quadrilateral meshes with and 
without midpoints in 2D and 3D. This work was supported by The Ministry of Education, Youth and 
Sports (of the Czech Republic) from the National Programme of Sustainability (NPU II), project 
“IT4Innovations excellence in science - LQ1602”. 


Appendix 

A. Simplified constitutive handling for the plane strain problem 

The results of Section and are, of course, valid also for the plane strain problem. Nevertheless, 
in this case, one can simplify the forms of eigenprojections and their derivatives since we work only 
on the subspace of containing trial tensors in the form 

( Vii Vi2 0 \ 
hi 2 V 22 0 . 

0 0 V33 / 

To distinguish the derivatives of functions dehned in Mps, we use the symbol V instead of V. Dehne 
the functions 


:= 2 

Vii + V22 + \/ivii - V22Y + 477^2 

^ 2 {'n) := ^ 

Vii + ri22- sj{rjii - V22f + 4?7^2 

^ 3 {'n) ■= V 33 



in Mpp. Then rji = i = 1,2,3, are the eigenvalues of 77 . These values are not ordered in 

general. We only know that > 7 ) 2 - Further, dehne 




hii 

V12 

M ~ 1 

f 1 

0 


Tll2 

V22 

0 , /:= 

0 

1 

0 

0 

0 

0 / 

1 0 

0 

0 / 




0 0 0 \ 

0 0 0 , 

0 0 1 / 


EM 


'n-ml 

vi-m ’ 


i, 


Vi > V2 
fji = m 


E2{'n) .= I - EM 
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and 


Ei(r/) 


[I — El ^ El — E2 ® E2 ], 

O, 


Vi > V 2 

vi = m 


^ 2 {'n) := -Ei(r/), E 3 (? 7 ) := O, 


where O denotes the zeroth fourth order tensor and [t\ijki = hJy k,l = 1, 2 , otherwise [t\ijki = 0. 

Clearly, Vu 3 {ri) = 'Du 3 {rj) = E 3 {r]). If fji > 7)2 then 

Ei{r]) = 'Dui{r]), Ei{r]) = 'bEi{r]) in Rps, i = 1,2. 


It is worth mentioning that these formulas need not hold in in general. Similar formulas are 
also introduced in [21 Appendix A], 

Now, it is necessary to reorder the eigenvalues of r/ G Rps- Denote the ordered eigenvalues as 
ViyV 2 ,V 3 ^ i- 6 -; hi •= max{? 7 i,? 73 } and rjs := min{r 72 ,h 3 }- Consequently, we reorder the functions a)*, 
Ei, Ej, i = 1,2,3, leading to the functions Ui, Ei, Ej, i = 1,2,3. To complete the notation, one can 
easily set 

Ei2{r)) := ^^1(77) + E2{r)), E 23 {v) ■= -£^2(^7) + -^3(77) V77 G Rps- 

Finally, one can straightforwardly use the functions uji, Ei, Ej, 7 = 1,2, 3, E 12 and E 23 within 
Section when the plane strain assumptions are considered. 


B. Algebraic representation of second and fourth order tensors 


Within our implementation, we use the standard algebraic representation of stress and strain sec¬ 
ond order tensors specihed below but a little bit different representation of fourth order tensors in 
comparison to [21 Appendix D]. We assume that a fourth order tensor C represents a linear mapping 
from into R^y^- Therefore, the components [C\ijki = Cijki of C satisfy 

CijkiVki = Y V77 G Rly^, riki = ['n]ki. 

k,l k,l 


The choice rjki = Smk^ni + ^nk^mi implies that 77 G R^y^ for any m, n = 1, 2, 3 and 


Cijmn T Cijnm C iimn + c 


-'tjnm 


jimn 


jinm 


'ii,i, m,n = 1, 2, 3. 


(B.l) 


Notice that in [21 Appendix D], the stronger assumptions on the components are required: Cijmn = 

r —c —c 

^ jimn 

We distinguish two cases: the 3D problem and its plane strain reduction. 


The 3 D problem 


Let T, 77 G R-lym denote stress and strain tensors, respectively. Then they are represented by vectors 
t = (rii,r 22 ,r 33 ,ri 2 ,r 23 ,ri 3 )'^ and n = ( 7711 , 7722 , 7733 , 27712 , 27723 , 27713 )^ where nj and 77 ^- are the com¬ 
ponents of T, and 77, respectively. Clearly, r : 77 = t ■ n. A fourth order tensor C is represented by 
matrix C G Since fourth order tensors are applied on strain tensors within the implementation, 

we require that 


77 : C : £ = n ■ Ce (B.2) 

holds for any strain tensors 77 and £. Here, n and e denote the algebraic counterparts of 77 and e, 
respectively. From (B.l) and (B.2), one can derive that 


C = 


/ 

Ciiii 

C1122 

Cl 133 

^ [Cl 112 

+ 

C1121] 

1 [C1123 

+ 

C1132] 

|[Ciii3 

+ 

C1131] 

\ 


C2211 

C2222 

C2233 

^[C2212 

+ 

C2221] 

|[C2223 

+ 

C2232] 

|[C2213 

+ 

C2231] 



C3311 

C3322 

C3333 

^ [C3312 

+ 

C3321] 

1 [C3323 

+ 

C3332] 

1 [C3313 

+ 

C3331] 



C1211 

C1222 

C1233 

^ [C1212 

+ 

C1221] 

|[Ci 223 

+ 

C1232] 

|[Ci213 

+ 

C1231] 



C2311 

C2322 

C2333 

^[C2312 

+ 

C2321] 

|[C2323 

+ 

C2332] 

|[C2313 

+ 

C2331] 


V 

C1311 

C1322 

C1333 

^ [C1312 

+ 

C1321] 

|[Ci323 

+ 

C1332] 

1 [Cl313 

+ 

C1331] 

/ 
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Indeed, the choices Ski = ^{ 62^1 + ^ 2 ^ 1 ), Vij = \{^uS 2 j + S 2 i 6 ij) imply e = (0, 0,0, 0,1, 0)^ and 
n = (0, 0, 0,1, 0, 0)^. Hence, 

[C]i^ = n ■ Ce 7 ] : C : e = -[6*1223 + 6*1232 + 6*2123 + 6*2132] 2 [^1223 + 6*1232] • 

Similarly, one can derive the forms of other components of C. Notice that the algebraic representation 
of C is more general than in [21 Appendix D], 

We introduce three examples useful for the Mohr-Coulomb model: 

1. Let C = I. Then Cijki = SikSji and C = diag(l, 1,1,1/2,1/2,1/2). Notice that the same matrix 

is derived in [21 Appendix D] although the tensor I 5 , = ^{SikSji + SuSjk), is used there 

instead of I. 

2. Let C = T ® cr where cr, r are arbitrary chosen stress tensors. Denote s and t as the algebraic 
counterparts to cr, r, respectively. Then C = st^. 

3. Let C = Then Cijki = Sikr]ij + 5jiT]ik and 


( 2 ? 7 ii 

0 

0 

V12 

0 

hl 3 \ 

0 

27722 

0 

V12 

V23 

0 

0 

0 

27733 

0 

V23 

Vl 3 

V12 

V12 

0 

l[Vll + V 22 ] 

hl 3 

\ri23 

0 

V23 

V23 

CO 

^[h 22 + ^733] 

\rii2 

\ Vl 3 

0 

Vl 3 

1^23 

lVl2 

\[jll+h 33 ] / 


The plane strain problem 

Let T and 77 denote stress and strain second order tensors, respectively. Then they are represented 
by the vectors t = (th, r22, t'i 2, rsa)^ and n = ( 7711 , 7722 ,2?7i2, 7733 )^ where Tij and rjij are components 
of T, and r], respectively. Clearly, r : 77 = t ■ n. Notice that the component 7733 vanishes for the 
strain tensor but not for the plastic strain tensor. 

The fourth order tensor C can be represented by matrix C G Similarly as for the 3D 

problem, one can derive that 


^ 6*1111 

6*1122 

1 [6*1112 + 6*1121] 

6*1133 

\ 

6*2211 

6*2222 

1 [6*2212 + 6*2221] 

6*2233 


6*1211 

6*1222 

1 [6*1212 + 6*1221] 

6*1233 


V 6*3311 

6*3322 

1 [ 6*3312 + 6*3321] 

6*3333 

/ 


Finally, it is worth mentioning that for assembling the tangent stiffness matrix, it is sufficient to save 
only the components (C')jj where i,j = 1,2,3. 


References 

[ 1 ] Sysala S, Cermak M, Koudelka T, Kruis J, Zeman J, Blaheta R. Subdifferential-based implicit 
return-mapping operators in computational plasticity. ZAMM-Journal of Applied Mathematics 
and Mechanics/Zeitschrift fiir Angewandte Mathematik und Mechanik 2016; . 

[2] de Souza Neto EA, Peric D, Owen DRJ. Computational Methods for Plasticity. Wiley-Blackwell, 
2008. 


24 






[3] Clausen J, Damkilde L, Andersen LV. Robust and efficient handling of yield surface discontinu¬ 
ities in elasto-plastic finite element calculations. Engineering Computations 2015; 32(6):1722- 
1752. 

[4] Lin C, Li YM. A return mapping algorithm for unihed strength theory model. International 
Journal for Numerical Methods in Engineering 2015; 104(8):749-766. 

[5] Larsson R, Runesson K. Implicit integration and consistent linearization for yield criteria of the 
mohr-coulomb type. Mechanics of Cohesive-frictional Materials 1996; l(4):367-383. 

[6] Borja RI, Sama KM, Sanz PF. On the numerical integration of three-invariant elastoplastic 
constitutive models. Computer Methods in Applied Mechanics and Engineering 2003; 192(9- 
10):1227-1258. 

[7] Karaoulanis FE. Implicit numerical integration of nonsmooth multisurface yield criteria in the 
principal stress space. Arch Computat Methods Eng 2013; 20(3):263-308. 

[8] Koiter WT. Stress-strain relations, uniqueness and variational theorems for elastic-plastic ma¬ 
terials with a singular yield surface. Quarterly of Applied Mathematics 1953; ll(3):350-354. 

[9] de Borst R. Integration of plasticity equations for singular yield functions. Computers & Struc¬ 
tures 1987; 26(5):823-829. 

[10] Han W, Reddy BD. Plasticity: mathematical theory and numerical analysis. Springer-Verlag, 
1999. 

[11] Berga A. Mathematical and numerical modeling of the non-associated plasticity of soils—part 1: 
The boundary value problem. International Journal of Non-Linear Mechanics 2012; 47(1):26- 
35. 

[12] Simo JC, Hughes TJ. Computational inelasticity. Springer Science & Business Media, 2006. 

[13] Starman B, Halilovic M, Vrh M, Stok B. Consistent tangent operator for cutting-plane algorithm 
of elasto-plasticity. Computer Methods in Applied Mechanics and Engineering 2014; 272:214- 
232. 

[14] Clausen J, Damkilde L, Andersen L. Efficient return algorithms for associated plasticity 
with multiple yield planes. International Journal for Numerical Methods in Engineering 2006; 
66(6):1036-1059. 

[15] Abbo A, Lyamin A, Sloan S, Hambleton J. A C2 continuous approximation to the 
Mohr-Coulomb yield surface. International Journal of Solids and Structures 2011; 48(21):3001- 
3010. 

[16] Borja RI. Plasticity. Springer, 2013. 

[17] Cermak M, Kozubek T, Sysala S, Valdman J. A TFETI domain decomposition solver for elasto¬ 
plastic problems. Applied Mathematics and Computation 2014; 231:634-653. 

[18] Gruber PC, Valdman J. Solution of one-time-step problems in elastoplasticity by a slant Newton 
method. SIAM J. Sci. Comput. 2009; 31 (2): 1558-1580. 

[19] Sauter M, Wieners C. On the superlinear convergence in computational elasto-plasticity. Com¬ 
puter Methods in Applied Mechanics and Engineering 2011; 200(49-52):3646-3658. 


25 



[20] Sysala S. Application of a modified semismooth Newton method to some elasto-plastic problems. 
Mathematics and Computers in Simulation 2012; 82(10):2004-2021. 

[21] Sysala S. Properties and simplihcations of constitutive time-discretized elastoplastic operators. 
ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift fiir Angewandte Mathe- 
matik und Mechanik 2014; 94(3):233-255. 

[22] Carlson DE, Roger A. The derivative of a tensor-valued function of a tensor. The Quarterly of 
Applied Mathematics 1986; 44:409-423. 

[23] de Borst R, Crisheld MA, Remmers JJC, Verhoosel CV. Non-Linear Finite Element Analysis 
of Solids and Structures. Wiley-Blackwell, 2012. 

[24] Chen WF, Liu X. Limit analysis in soil mechanics. Elsevier, 2012. 

[25] Sysala S, Haslinger J, Hlavacek I, Cermak M. Discretization and numerical realization of contact 
problems for elastic-perfectly plastic bodies. PART I - discretization, limit analysis. ZAMM- 
Journal of Applied Mathematics and Mechanics/Zeitschrift fiir Angewandte Mathematik und 
Mechanik 2015; 95(4):333-353. 

[26] Cermak M, Haslinger J, Kozubek T, Sysala S. Discretization and numerical realization of con¬ 
tact problems for elastic-perfectly plastic bodies. PART II - numerical realization, limit analysis. 
ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift fiir Angewandte Mathe¬ 
matik und Mechanik 2015; 95(12):1348-1371. 

[27] Haslinger J, Repin S, Sysala S. A reliable incremental method of computing the limit load 
in deformation plasticity based on compliance: Continuous and discrete setting. Journal of 
Computational and Applied Mathematics 2016; 303:156-170. 

[28] Haslinger J, Repin S, Sysala S. Guaranteed and computable bounds of the limit load for vari¬ 
ational problems with linear growth energy functionals. Applications of Mathematics 2016; 
61(5):527-564. 

[29] Ruszczyhski AP. Nonlinear optimization. Princeton university press, 2006. 

[30] Sysala S, Cermak M. Experimental matlab code for the slope stability benchmark - SS-MC-NP- 
3D, SS-MC-NH, SS-MC-NP-Acontrol 2016. URL www.ugn. cas . cz/?p=publish/output .php, 
(or www.ugn.cas.cz - Publications - Other outputs - SS-MC-NP-3D, SS-MC-NH, SS-MC-NP- 
Acontrol) . 


26 


